index

1. Data

Show the code

# Library
library(dplyr)
library(ggplot2)

# load data 
## Isotope data 
isotope_data <-  utils::read.csv(here::here("data", "isotopic_data_2021.csv"), sep = ";", header = T,dec = ",")

isotope_data_fish <- isotope_data %>%
  # only fish 
  filter(taxon == "Fish")%>%
  mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))%>%
  arrange(species)

## trawling data 
trawling_data_evhoe21 <-  utils::read.csv(here::here("data", "trawling_data_evhoe_2021.csv"), sep = ";", header = T, dec = ".")

# density distribution from biomass value 
density_distribution <- trawling_data_evhoe21%>%
  select(Nom_Scientifique, Tot_V_HV, Code_Station)%>%
  #selection of mesopelagic trawling 
  filter(Code_Station%in% c("Z0524", "Z0518", "Z0512", "Z0508", 
                            "Z0503", "Z0497", "Z0492"))%>%
  #selection of species sampled for isotope 
  filter(Nom_Scientifique%in%c("Arctozenus risso", "Argyropelecus olfersii", "Benthosema glaciale",
                               "Cyclothone", "Lampanyctus crocodilus", "Lampanyctus macdonaldi",
                               "Lestidiops sphyrenoides", "Maulisia argipalla", "Maurolicus muelleri",
                               "Melanostigma atlanticum", "Myctophum punctatum", "Notoscopelus bolini",
                               "Notoscopelus kroyeri", "Searsia koefoedi", "Serrivomer beanii",
                               "Xenodermichthys copei"))%>%
  mutate(Nom_Scientifique=recode(Nom_Scientifique, "Cyclothone"="Cyclothone spp."))%>%
  mutate(trawling_depth= case_when(Code_Station %in% c("Z0508") ~25,
                                   Code_Station %in% c("Z0492") ~370,
                                   Code_Station%in% c("Z0512") ~555,
                                   Code_Station %in% c("Z0503") ~715,
                                   Code_Station %in% c("Z0518") ~1000,
                                   Code_Station %in% c("Z0524") ~1010,
                                   Code_Station %in% c("Z0497") ~1335))%>%
  distinct()%>%
  group_by(Nom_Scientifique)%>%
  mutate(sum_sp=sum(Tot_V_HV))%>%
  ungroup()%>%
  group_by(trawling_depth, Nom_Scientifique)%>%
  mutate(pourcentage_bio=sum(Tot_V_HV/sum_sp*100))%>%
  # Selection of trawling depth 
  select(Nom_Scientifique, trawling_depth, pourcentage_bio)%>%
  # to have a round number to be able to multiply it afterwards 
  mutate(across(pourcentage_bio, round, 0)) %>%
  mutate(n_bio = as.integer(pourcentage_bio))%>%
  select(-c(pourcentage_bio))%>%
  tidyr::uncount(n_bio)

2. Model deterministic or stochastic ?

Objectif: déterminer si les espèces cooccurrentes au sein d’une couche de profondeur présentent des niches trophiques compatibles avec un modèle déterministe d’aménagement de niche ou un modèle stochastique (dû au hasard).

Hypothèses:

• Si modèle déterministe: la composition et la proportion des espèces observées dans un assemblage n’est pas dû au hasard, la concurrence pour les ressources pousse les espèces à devenir moins semblables dans leur utilisation des ressources, ce qui réduira la taille de leur niche et favorisera donc leur coexistence en limitant leurs similitudes

• Si modèle stochastique: composition et l’abondance des espèces d’un assemblage est dû au hasard, la coexistence se produit par le biais de l’équivalence écologique (les espèces possèdent des caractéristiques écologiques similaires et n’ont pas d’avantage concurentiel les unes sur les autres)

→ Pour cela :

• modèle null (stochastique) = distribution attendue lorsque les modèles résultent de processus stochastiques (en rééchantillonnant les données isotopiques), donc un assemblage organisé par des processus stochastiques devrait présenter des niches trophiques avec des positions similaires dans l’espace isotopique, une grande taille de niche et des niveaux élevés de chevauchement entre les espèces.

• Si différence avec le modèle null = assemblage oragnisé par des intéractions déterministes entre espèces de compétition et de prédation. Donc position des niches différenciée entre les espèces, taille des niches relativement petite par rapport à tout l’espace isotopique occupé de l’assemblage, et chevauchement relativement faible entre les espèces.

Calcul de différents indices:

• Taille ellipse à 40%
  
• Mesures des chevauchements interspécifiques

Taille des niches

• ellipses à 40%

niches_area_model

conclusion : modèle déterministe (niches observées plus petites que si distribution des espèces au hasard) pour la plupart des couches de profondeur, sauf couche épipélagique (et station près du fond). Ressources non limitantes dans la couche épipélagique (et près du fond) donc les espèces ne sont pas forcées de se ségréger ?

Chevauchements

• Calcul de la matrix de chevauchement entre les espèces (ellipses à 40%), puis somme de tous les chauvauchements.

conclusion: - les chevauchements observés sont significativement différents que si modèle de distribution au hasard, modèle déterministe. (mais toutes les espèces présentes dans un assemblage non échantillonées, pb?) - Lower-mesopelagic: nombre d’espèces échantillonnées le plus important (n=9) et valeur observée la plus petite par rapport à la distribution des valeurs permutées: pour le nombre d’espèces présentes peu d’overlap? (et inversement pour bottom proximity?)

Indices de diversité isotpiques

• sans prise en compte de la biomasse

conclusion: modèle stochastique. Sans pondération de la biomasse toutes les espèces ont le même poids dans l’analyse.

A garder dans le papier? :

• Taille des ellipses
• chevauchements
• range range carbone et azote

Comme modèle déterministe, quelles sont les forces évolutives/les axes de ségragtion qui rentrent en jeu pour expliquer cette structure?

• Calcul des indices de disversité isotopiques pondéré par la biomasse ? La structure de la communauté peut être aussi expliquée par la répartition de la biomassse des espèces au sein de l’espace isotopique ?
• Autre axe de ségrégation le long du gradient de niveau trophique
• Le long du gradient de profondeur (analyse cluster guildes trophiques avec mise en lien de la distribution en profondeur des espèces)

Show the code

isotope_data_sum <- isotope_data%>%
  group_by(species)%>%
  mutate(mean_d13c=mean(d13c),
         sd_d13c=sd(d13c),
         mean_d15n=mean(d15n),
         sd_d15n=sd(d15n),
         n=n())%>%
  select(species, mean_d13c, sd_d13c, mean_d15n, sd_d15n, n, number_ind_pool, family)%>%
  distinct()%>%
  mutate(across(where(is.numeric), round, 2))

htmltools::tagList(DT::datatable(isotope_data_sum))

3. Isotopic niches

Ellipses

• ellipses at 40%
• Δ \(\delta\)¹³C = 2.36‰
• Δ \(\delta\)¹⁵N = 5.94‰

Show the code

niche_plot_community <- isotope_data%>%
  mutate(species= gsub("_"," ", species))%>%
    mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# colors selection 
nichecol_krill <- c("#E4A33A", "#F67451", "#D664BE", "#3DA5D9", "#94b3ae",
                    "#18206F", "#FD151B", "#049A8F", "#072AC8", "black", "purple", 
                    "#d193f7", "#d8c2ab", "#678FCB", "#A63A49", "#00547A", "#6B54A0")

# size of the ellipse 
p.ell<-0.40
 
# plot
ggplot(data = niche_plot_community, 
                         aes(x = d13c, 
                             y = d15n)) + 
  geom_point(aes(color = species, shape= taxon), size = 2) +
  scale_color_manual(values=nichecol_krill)+  
  scale_fill_manual(values=nichecol_krill)+
  scale_shape_manual(values= c(19, 3))+
  scale_x_continuous(expression({delta}^13*C~'\u2030'), limits = c(-21, -18.5)) +
  scale_y_continuous(expression({delta}^15*N~'\u2030'), limits = c(7, 14))+
  stat_ellipse(aes(group = species, fill = species, color = species), 
               alpha = 0.2, level = p.ell,linewidth = 0.7, type = "norm", geom = "polygon")+
  theme_light()+
  theme(legend.text = element_text(size=15),
        legend.title = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15))+
  labs(shape="Taxon", col= "Species", fill="Species")

Show the code

#ggsave("niches_community.png", path = "figures", dpi = 700, height = 10, width = 12)

Niche areas

Show the code

# Model initialization ----
# options for running jags
parms <- list()
parms$n.iter <- 2 * 10^4   # number of iterations to run the model for
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10     # thin the posterior by this many
parms$n.chains <- 2        # run this many chains
parms$save.output = FALSE
parms$save.dir = tempdir()

# define the priors
priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3

# Niche area ----
community_siber <- isotope_data%>%
  select(-group)%>%
  mutate(group = species, 
         # All species in the same community (no segregation by depth)
         community = "",
         iso1=d13c,
         iso2 =d15n,
         .keep = "unused")%>%
  select(group, community, iso1, iso2)%>%
  relocate(iso1, .before = group)%>%
  relocate(iso2, .after = iso1)%>%
  # Reorder species by SEA
  mutate(group_order = forcats::fct_relevel(group, c("Searsia_koefoedi", "Lampanyctus_crocodilus",
                                              "Myctophum_punctatum", "Benthosema_glaciale",
                                              "Xenodermichthys_copei", "Lampanyctus_macdonaldi",
                                              "Serrivomer_beanii", "Cyclothone spp.", "Argyropelecus_olfersii",
                                              "Melanostigma_atlanticum", "Arctozenus_risso", "Maulisia_argipalla",
                                              "Notoscopelus_kroyeri", "Lestidiops_sphyrenoides",
                                              "Notoscopelus_bolini", "Maurolicus_muelleri", "Meganyctiphanes_norvegica"))) %>%
  arrange(group_order)%>%
  select(-group)%>%
  rename("group"="group_order")%>%
  relocate(community, .after = group)

# create the siber object
community_siber_obj <- SIBER::createSiberObject(community_siber)

# Calculate summary statistics for each group: TA, SEA and SEAc
group_ML <- SIBER::groupMetricsML(community_siber_obj)

# fit the ellipses which uses an Inverse Wishart prior
# on the covariance matrix Sigma, and a vague normal prior on the 
# means. Fitting is via the JAGS method.
ellipses_posterior_community <- SIBER::siberMVN(community_siber_obj, parms, priors)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 14
   Unobserved stochastic nodes: 3
   Total graph size: 29

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 120
   Unobserved stochastic nodes: 3
   Total graph size: 135

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 57
   Unobserved stochastic nodes: 3
   Total graph size: 72

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 97
   Unobserved stochastic nodes: 3
   Total graph size: 112

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 26
   Unobserved stochastic nodes: 3
   Total graph size: 41

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 41
   Unobserved stochastic nodes: 3
   Total graph size: 56

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 43
   Unobserved stochastic nodes: 3
   Total graph size: 58

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 14
   Unobserved stochastic nodes: 3
   Total graph size: 29

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 60
   Unobserved stochastic nodes: 3
   Total graph size: 75

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 12
   Unobserved stochastic nodes: 3
   Total graph size: 27

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 35
   Unobserved stochastic nodes: 3
   Total graph size: 50

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Show the code

SEA_B_community <- SIBER::siberEllipses(ellipses_posterior_community)

# colors 
clrs_com <- matrix(c( "#e4c3c8", "#c98891", "#A63A49",
                      "#dee8e6", "#bed1ce", "#94b3ae", 
                      "#F5E9FD", "#E7C9FB", "#d193f7",
                      "#f2d0eb", "#e6a2d8", "#D664BE",
                      "#d2cbe2", "#a698c6", "#6B54A0",
                      "#b9bcd3", "#5d629a", "#18206F",
                      "#b2cbd7", "#6698af", "#00547A",
                      "#c4e4f3", "#8ac9e8", "#3DA5D9",
                      "#fcd5ca", "#f9ab96", "#F67451",
                      "#ECD7F9", "#D09CF1", "purple",
                      "#F9ECD7", "#F1D09C", "#E4A33A",
                      "#cceae8", "#68c2bb", "#049A8F",
                      "#c2d2ea", "#a3bbdf", "#678FCB",
                      "#fed0d1", "#fe8a8d", "#FD151B",
                      "#f3ece5", "#e7dacc", "#d8c2ab",
                      "#CDD4F4", "#8394E3", "#072AC8"), nrow = 3, ncol = 16)

# save plot in high quality 
#tiff("SEA.tiff", units="in", width=10, height=7, res=700)

# plot
SIBER::siberDensityPlot(SEA_B_community, xticklabels = colnames(group_ML),
                        ylab = expression("Standard Ellipse Area" ('\u2030' ^2) ),
                        bty = "L",
                        las = 1,
                        xlab = "",
                        #clr = clrs_com,
                        ylims = c(0,1.5))

# Add x's for the ML estimated SEA-c
points(1:ncol(SEA_B_community), group_ML[3,], col="red", pch = "x", lwd = 2)

Show the code

#dev.off()

cr.p <- c(0.95, 0.99) # vector of quantiles

# do similar to get the modes, taking care to pick up multimodal posterior
# distributions if present
SEA_B_community_modes <- lapply(
  as.data.frame(SEA_B_community), 
  function(x,...){tmp<-hdrcde::hdr(x)$mode},
  prob = cr.p, all.modes=T)

Overlaps

Show the code

# Prepare data 
mat_overlap_data <- isotope_data%>%
  filter(taxon=="Fish")%>%
  mutate(species= gsub("_"," ", species))%>%
  mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# Arrange species by their d15n values 
mat_overlap <- mat_overlap_data%>%
  group_by(species)%>%
  mutate(mean_d15n=mean(d15n))%>%
  arrange(mean_d15n)%>%
  as.data.frame()

# ellipse at 40%
test.elp_community <- rKIN::estEllipse(data= mat_overlap,  x="d13c", y="d15n", group="species", levels=40, smallSamp = TRUE)

[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"

Show the code

# Extract the area of each polygon
elp.area_community <- rKIN::getArea(test.elp_community)

# determine polygon overlap for all polygons
elp.olp_community <- rKIN::calcOverlap(test.elp_community)%>%
  #renames rownames
  mutate(OverlapID= gsub("_"," ", OverlapID))%>%
  mutate(OverlapID= gsub("40","", OverlapID))%>%
  mutate(across(where(is.numeric), round, 2))%>%
  tibble::column_to_rownames("OverlapID")%>%
  as.data.frame()

#renames colnames
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "_", replacement = " ")
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "40", replacement = "")

#plot 
ggcorrplot::ggcorrplot(elp.olp_community, lab = T, outline.color = "white", lab_size = 3, tl.cex = 10)+
  scale_fill_gradient2(limit = c(0,1), low = "white", high = "grey50", mid = "grey80", midpoint = 0.5)+
  labs(fill="Overlap value")+
  theme(axis.text = element_text(face="italic", size = 10),
        legend.text = element_text(size=10),
        legend.title = element_text(size=10),
        plot.background = element_rect(colour = "white"))

Show the code

#ggsave("matrix_overlap.png", path = "figures", dpi = 700, width = 8, height = 6)

4. Depth segregation

Cluster

• input data for cluster = overlap matrix
• depth distribution = from complete 2021 trawling data (not only individuals sampled for isotope)

Show the code

# calculate median depth by species with trawling data set 2021
mean_depth_sp <- density_distribution%>%
  group_by(Nom_Scientifique)%>%
  mutate(median_depth=median(trawling_depth))%>%
  select(Nom_Scientifique, median_depth)%>%
  distinct()%>%
  ungroup()%>%
  arrange(Nom_Scientifique)%>%
  rename(species=Nom_Scientifique)

# calcul overlap ellipse 40%
# test.elp_community <- rKIN::estEllipse(data= isotope_data_fish,  x="d13c", y="d15n", group="species", levels=40, smallSamp = TRUE)
# elp.olp_community <- rKIN::calcOverlap(test.elp_community)
# 
# df <- elp.olp_community%>%
#   tibble::column_to_rownames(var="OverlapID")%>%
#   as.data.frame()

# Gap statistic
factoextra::fviz_nbclust(elp.olp_community , kmeans, nstart = 25,  method = "gap_stat", nboot = 100, verbose = FALSE)+
  labs(subtitle = "Gap statistic method")

Show the code

res.km <- kmeans(scale(elp.olp_community ),5, nstart = 25)

# Réduction de dimension en utilisant l'ACP
res.pca <- prcomp(elp.olp_community [, -5],  scale = TRUE)
# Coordonnées des individus
ind.coord <- as.data.frame(factoextra::get_pca_ind(res.pca)$coord)
# Ajouter les clusters obtenus à l'aide de l'algorithme k-means
ind.coord$cluster <- factor(res.km$cluster)
# Percentage of variance explained by dimensions
eigenvalue <- round(factoextra::get_eigenvalue(res.pca), 1)
variance.percent <- eigenvalue$variance.percent

# Add species names and median depth 
ind.coord_depth <- cbind(ind.coord, mean_depth_sp)
ind.coord_depth$median_depth <- as.factor(ind.coord_depth$median_depth)

Cluster and depth distribution

Show the code

# dendrogram ----
dend <- elp.olp_community %>%
  dist() %>%
  hclust() %>%
  as.dendrogram()

#png("figures/dendrogram.png", units="in", width=6, height=4, res=700)
par(mar = c(1, 1, 1, 10))
dend %>%
  dendextend::set("labels_col",
      value = c("#86BBBD", "#ECA72C", "#4D85A8", "#9BABE8", "#D35D4A"),
      k = 5) %>%
  dendextend::set("branches_k_color",
      value = c("#86BBBD", "#ECA72C", "#4D85A8", "#9BABE8", "#D35D4A"),
      k = 5) %>%
  dendextend::set("labels_cex", 0.6) %>%
  dendextend::set("branches_lty", 2) %>%
  dendextend::set("leaves_bg") %>%
  dendextend::set("leaves_pch", 19) %>%
  dendextend::set("leaves_col",
    c(
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#ECA72C",
      "#ECA72C",
      "#ECA72C",
      "#ECA72C",
      "#4D85A8",
      "#4D85A8",
      "#4D85A8",
      "#9BABE8",
      "#9BABE8",
      "#D35D4A",
      "#D35D4A"
    )
  ) %>%
  plot(horiz = TRUE, axes = FALSE)

Show the code

#dev.off()

# Density plot----
# assign each species to a cluster 
density_distribution_cluster <- density_distribution%>%
  mutate(cluster=case_when(Nom_Scientifique%in% c("Argyropelecus olfersii", "Lampanyctus crocodilus",
                                                  "Benthosema glaciale","Myctophum punctatum","Serrivomer beanii")~4,
                           Nom_Scientifique%in% c("Lampanyctus macdonaldi", "Maulisia argipalla",
                                                  "Searsia koefoedi")~5,
                           Nom_Scientifique%in% c("Cyclothone spp.","Notoscopelus bolini", "Notoscopelus kroyeri",
                                                  "Melanostigma atlanticum")~3,
                           Nom_Scientifique%in% c("Xenodermichthys copei", "Maurolicus muelleri")~1,
                           Nom_Scientifique%in% c("Arctozenus risso","Lestidiops sphyrenoides")~2))%>%
  group_by(Nom_Scientifique)%>%
  arrange(desc(trawling_depth))

# Order in function of median depth
density_distribution_cluster$Nom_Scientifique = with(density_distribution_cluster, reorder(Nom_Scientifique, cluster, max))  

ggplot(density_distribution_cluster,
       aes(x = trawling_depth, y = Nom_Scientifique, group = Nom_Scientifique, 
           col=factor(cluster), fill=factor(cluster)))+ 
  scale_fill_manual(values = c("#D35D4A", "#9BABE8", "#ECA72C", "#86BBBD","#4D85A8"))+
  scale_color_manual(values = c("#D35D4A", "#9BABE8", "#ECA72C", "#86BBBD","#4D85A8"))+
  ggridges::stat_density_ridges(quantile_lines = TRUE, quantiles = 0.5 , alpha=0.4, size=0.7,
                                rel_min_height = 0.002, scale=1.2)+
  theme_light()+
  scale_y_discrete(position = "left")+
  scale_x_reverse(limits = c( 1400,0))+
  coord_flip()+
  ylab(label = "")+ xlab("Depth (m)")+
  theme(axis.text.y = element_text(size=10.5),
        axis.text.x = element_text(face="italic", size=10.5, angle=80, vjust = 0.5, hjust=0.5),
        axis.title.x = element_text(size=10.5),
        axis.title.y = element_text(size=10.5))+
  guides(fill="none", col="none", alpha="none")

Show the code

ggsave("density_plot.png", path = "figures", dpi = 600, height = 8, width = 10)

5. Isotopic diversity index

Definitions

• Isotopic divergence : tends to 1 when all points (or their weights) are located at the edge of the convex hull, i.e. when the oragnisms with the most extreme isotopic values dominate the food web

• Isotopic dispersion: equal to 1 when most points (or their weights) are far from the center of gravity of the point group, i.e. when organisms tend to have contrasting isotopic values

• Isotopic evenness: tends to 1 when all organisms are equally distributed in isotopic space

• Isotopic uniqueness: tends to 1 when most organisms (or their weights) are isolated in isotopic space, i.e. when most organisms (or those with the highest abundance) are isolated in isotopic space, their isotopic values are very different from the rest of the organisms

Representativeness of sampling

• Percentage of species biomass in each station

Show the code

# Species % biomass and abundance by station 
species_abundance <- trawling_data_evhoe21%>%
  # selection of mesopelagic trawl 
  filter(Code_Station%in% c("Z0524", "Z0518", "Z0512", "Z0508", 
                            "Z0503", "Z0497", "Z0492"))%>%
  #deletion when only genus name and genus already sampled in isotopy + A. carbo
  filter(!Nom_Scientifique%in%c("Cyclothone braueri","Cyclothone microdon", 
                                "Myctophidae", "Aphanopus carbo"))%>%
    mutate(
    depth_layer = case_when(
      Code_Station %in% c("Z0508") ~ "epipelagic",
      Code_Station %in% c("Z0492", "Z0512") ~ "upper-mesopelagic",
      Code_Station %in% c("Z0503", "Z0518") ~ "lower-mesopelagic",
      Code_Station %in% c("Z0497") ~ "bathypelagic",
      Code_Station %in% c("Z0524") ~ "bottom-proximity")) %>%
  select(depth_layer, Nbr, Nom_Scientifique, Tot_V_HV)%>%
  group_by(Nom_Scientifique, depth_layer)%>%
  mutate(nb_ind= sum(Nbr))%>%
  select(-Nbr)%>%
  distinct()%>%
  group_by(depth_layer)%>%
  mutate(sum_biomass_station=sum(Tot_V_HV))%>%
  ungroup()%>%
  group_by(depth_layer, Nom_Scientifique)%>%
  mutate(pourcentage_biomass_sp= Tot_V_HV/sum_biomass_station*100)%>%
  select(Nom_Scientifique, pourcentage_biomass_sp, nb_ind, depth_layer)%>%
  mutate(across(where(is.numeric), round, 2))
  
htmltools::tagList(DT::datatable(species_abundance))

Calculation

• Formatting of data and calculation of relative biomass within each station

Show the code

# Load data 
isotope_data_fish <- isotope_data %>% filter (species != "Meganyctiphanes_norvegica")

# sourcing the R functions from 'si_div' R script
#source("R/si_div.R")
source("R/si_div.R")

# Format indiviudal_si
individuals_si <- isotope_data_fish %>%
  select(individual_code, station, d13c, d15n, species_code) %>%
  rename(indiv_ID=individual_code,
         d13C= d13c,
         d15N =d15n,
         Species_code= species_code)%>%
  arrange(Species_code)

# species_status_biomass ----
# with all species sampled (not only species sampled for isotope)

# load total trawling data evhoe 2021
trawling_data_evhoe21 <-  utils::read.csv(here::here("data", "trawling_data_evhoe_2021.csv"), sep = ";", header = T, dec = ".")


species_status_biomass <- trawling_data_evhoe21 %>%
  # selection of mesopelagic trawl
  filter(Code_Station %in% c("Z0524", "Z0518", "Z0512", "Z0508",
                             "Z0503", "Z0497", "Z0492")) %>%
  #deletion when only genus name and genus already sampled in isotope + A. carbo
  filter(
    !Nom_Scientifique %in% c(
      "Cyclothone braueri",
      "Cyclothone microdon",
      "Myctophidae",
      "Aphanopus carbo",
      "Lampanyctus"
    )
  ) %>%
  # group by depth layer 
  mutate(
    Status = case_when(
      Code_Station %in% c("Z0508") ~ "epipelagic",
      Code_Station %in% c("Z0492", "Z0512") ~ "upper-mesopelagic",
      Code_Station %in% c("Z0503", "Z0518") ~ "lower-mesopelagic",
      Code_Station %in% c("Z0497") ~ "bathypelagic",
      Code_Station %in% c("Z0524") ~ "bottom-proximity"
    )
  ) %>%
  select(Status,
         Tot_V_HV,
         Nom_Scientifique,
         Code_Espece_Campagne) %>%
  distinct() %>%
  # sum species biomass by depth layer
  group_by(Nom_Scientifique, Status) %>%
  mutate(biomass_sp = sum(Tot_V_HV)) %>%
  select(-Tot_V_HV) %>%
  distinct() %>%
  # sum of total biomass by depth layer
  group_by(Status) %>%
  mutate(biomass_tot = sum(biomass_sp)) %>%
  # relative biomass of each species by station
  mutate(rel_biomass = biomass_sp / biomass_tot * 100) %>%
  select(-c(biomass_sp, biomass_tot)) %>%
  # selection of species sampled for isotopy in each depth
  filter(
    Status == "epipelagic" &
      Nom_Scientifique %in% c(
        "Lestidiops sphyrenoides",
        "Maurolicus muelleri",
        "Myctophum punctatum",
        "Notoscopelus kroyeri",
        "Notoscopelus bolini"
      ) |
      Status == "upper-mesopelagic" &
      Nom_Scientifique %in% c(
        "Arctozenus risso",
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Myctophum punctatum",
        "Notoscopelus kroyeri",
        "Xenodermichthys copei"
      ) |
      Status == "lower-mesopelagic" &
      Nom_Scientifique %in% c(
        "Arctozenus risso",
        "Argyropelecus olfersii",
        "Cyclothone",
        "Benthosema glaciale",
        "Lampanyctus crocodilus",
        "Maulisia argipalla",
        "Searsia koefoedi",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      ) |
      Status == "bathypelagic" &
      Nom_Scientifique %in% c(
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Lampanyctus macdonaldi",
        "Myctophum punctatum",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      ) |
      Status == "bottom-proximity" &
      Nom_Scientifique %in% c(
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Melanostigma atlanticum",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      )
  ) %>%
  mutate(Species_code = tolower(Code_Espece_Campagne)) %>%
  rename(Species_name = Nom_Scientifique) %>%
  select(-Code_Espece_Campagne) %>%
  relocate(Status, .after = rel_biomass) %>%
  relocate(Species_code, .after = Species_name) %>%
  arrange(Species_name) %>%
  distinct()

htmltools::tagList(DT::datatable(species_status_biomass))

• percentage of biomass sampled for isotopy in each depth layer

Show the code

biomass_sampled <- species_status_biomass%>%
  group_by(Status)%>%
  mutate(sum_biomass = sum(rel_biomass))%>%
  select(sum_biomass, Status)%>%
  distinct()

htmltools::tagList(DT::datatable(biomass_sampled))

Epipelagic

_ 25m (Z0508)

Show the code

# 25m Z0508 ----
individuals_si_epi <- individuals_si%>%
  filter(station =="Z0508")%>%
  select(-station)

status_biomass_epi <- species_status_biomass%>%
  filter(Status=="epipelagic")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_epi<-data.frame(group=individuals_si_epi[,"Species_code"], individuals_si_epi)
mean_si_species_epi<-meanSI_group(individuals_si_epi)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_epi[,"sd_d13C"]/mean_si_species_epi[,"d13C"], CV_d15N=mean_si_species_epi[,"sd_d15N"]/mean_si_species_epi[,"d15N"] )

               CV_d13C    CV_d15N
lesti-sph -0.012948117 0.03266115
maur-mue  -0.005224151 0.05249590
myct-pun  -0.009814716 0.04047466
noto-bol  -0.009178283 0.02755544
noto-kro  -0.012155251 0.02163464

Show the code

# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_epi)==status_biomass_epi[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_epi <-data.frame(mean_si_species_epi[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_epi[,"rel_biomass"], Status=status_biomass_epi[,"Status"], latin_name=status_biomass_epi[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_epi<-scaleSI_range01(data_fish_epi)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_epi<-IDiversity(cons=data_fish_scl_epi, weight=data_fish_scl_epi[,c("rel_biomass")], nm_plot="1_epipelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_epi,3)) 

epipelagic_d13C_d15N

Upper mesopelagic

• 370m (Z0492) & 555m (Z0512)

Show the code

individuals_si_upm <- individuals_si%>%
  filter(station%in%c("Z0492", "Z0512"))%>%
  select(-station)

status_biomass_upm <- species_status_biomass%>%
  filter(Status=="upper-mesopelagic")

# computing mean Stable Isotope values for each species
individuals_si_upm<-data.frame(group=individuals_si_upm[,"Species_code"], individuals_si_upm)
mean_si_species_upm<-meanSI_group(individuals_si_upm)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_upm[,"sd_d13C"]/mean_si_species_upm[,"d13C"], CV_d15N=mean_si_species_upm[,"sd_d15N"]/mean_si_species_upm[,"d15N"] )

              CV_d13C    CV_d15N
arct-ris -0.010421638 0.03438256
argy-olf -0.009906681 0.04283576
lamp-cro -0.015534945 0.07982239
myct-pun -0.019927081 0.04074265
noto-kro -0.013115883 0.02239182
xeno-cop -0.012364479 0.07722547

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_upm <-data.frame(mean_si_species_upm[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_upm[,"rel_biomass"], Status=status_biomass_upm[,"Status"], latin_name=status_biomass_upm[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_upm<-scaleSI_range01(data_fish_upm)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_upm<-IDiversity(cons=data_fish_scl_upm, weight=data_fish_scl_upm[,c("rel_biomass")], nm_plot="2_upper_meso")

# printing results
result <- as.data.frame(round(ID_scl_ab_upm,3)) 

370m

Lower mesopelagic

• 715m (Z0503)
• 1000m (Z0518)

Show the code

individuals_si_lwm <- individuals_si%>%
  filter(station%in%c("Z0503", "Z0518"))%>%
  select(-station)

status_biomass_lwm<- species_status_biomass%>%
  filter(Status=="lower-mesopelagic")

# computing mean Stable Isotope values for each species
individuals_si_lwm<-data.frame(group=individuals_si_lwm[,"Species_code"], individuals_si_lwm)
mean_si_species_lwm<-meanSI_group(individuals_si_lwm)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_lwm[,"sd_d13C"]/mean_si_species_lwm[,"d13C"], CV_d15N=mean_si_species_lwm[,"sd_d15N"]/mean_si_species_lwm[,"d15N"] )

              CV_d13C    CV_d15N
arct-ris -0.013383792 0.02784671
argy-olf -0.010245107 0.03092443
bent-gla -0.015512985 0.06444602
cycl-otz -0.009264821 0.04946692
lamp-cro -0.020344506 0.04819340
maul-arg -0.009604780 0.03191368
sear-koe -0.022852986 0.05388211
serr-bea -0.012134196 0.06183851
xeno-cop -0.012017717 0.06196185

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_lwm <-data.frame(mean_si_species_lwm[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_lwm[,"rel_biomass"], Status=status_biomass_lwm[,"Status"], latin_name=status_biomass_lwm[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_lwm<-scaleSI_range01(data_fish_lwm)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_lwm<-IDiversity(cons=data_fish_scl_lwm, weight=data_fish_scl_lwm[,c("rel_biomass")], nm_plot="3_lower-mesopelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_lwm,3)) 

715m

Bathypelagic

• 1335m (Z0497)

Show the code

individuals_si_bathy <- individuals_si%>%
  filter(station=="Z0497")%>%
  select(-station)

status_biomass_bathy <- species_status_biomass%>%
  filter(Status=="bathypelagic")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_bathy<-data.frame(group=individuals_si_bathy[,"Species_code"], individuals_si_bathy)
mean_si_species_bathy<-meanSI_group(individuals_si_bathy)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_bathy[,"sd_d13C"]/mean_si_species_bathy[,"d13C"], CV_d15N=mean_si_species_bathy[,"sd_d15N"]/mean_si_species_bathy[,"d15N"] )

              CV_d13C    CV_d15N
argy-olf -0.009889227 0.01043410
lamp-cro -0.018677334 0.05649516
lamp-mac -0.022058396 0.02734756
myct-pun -0.024242014 0.04136229
serr-bea -0.014132768 0.05883833
xeno-cop -0.012585074 0.05902294

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_bathy <-data.frame(mean_si_species_bathy[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_bathy[,"rel_biomass"], Status=status_biomass_bathy[,"Status"], latin_name=status_biomass_bathy[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_bathy<-scaleSI_range01(data_fish_bathy)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_bathy<-IDiversity(cons=data_fish_scl_bathy, weight=data_fish_scl_bathy[,c("rel_biomass")], nm_plot="4_bathypelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_bathy,3)) 

bathypelagic

Near bottom

• Z0524

Show the code

# 1010 Z0524 ----
individuals_si_nb <- individuals_si%>%
  filter(station=="Z0524")%>%
  select(-station)

status_biomass_nb <- species_status_biomass%>%
  filter(Status=="bottom-proximity")

# computing mean Stable Isotope values for each species
individuals_si_nb<-data.frame(group=individuals_si_nb[,"Species_code"], individuals_si_nb)
mean_si_species_nb<-meanSI_group(individuals_si_nb)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_nb[,"sd_d13C"]/mean_si_species_nb[,"d13C"], CV_d15N=mean_si_species_nb[,"sd_d15N"]/mean_si_species_nb[,"d15N"] )

              CV_d13C    CV_d15N
argy-olf -0.016533711 0.03529939
lamp-cro -0.015886548 0.04471856
mela-atl -0.009854421 0.04078758
serr-bea -0.012700557 0.05654679
xeno-cop -0.009943149 0.04105076

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_nb <-data.frame(mean_si_species_nb[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_nb[,"rel_biomass"], Status=status_biomass_nb[,"Status"], latin_name=status_biomass_nb[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_nb<-scaleSI_range01(data_fish_nb)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_nb<-IDiversity(cons=data_fish_scl_nb, weight=data_fish_scl_nb[,c("rel_biomass")], nm_plot="5_near_bottom")

# printing results
result <- as.data.frame(round(ID_scl_ab_nb,3)) 

bathypelagic

Spider web

Show the code

# load library
library(tidyr)

# crate a data frame with all index value
trophic_indices <- data.frame(
  group = c("25m", "370m", "555m", "715m", "1000m", "1335m", "1010m-bottom proximity"),
  IDiv = c(0.952, 0.918, 0.845, 0.928, 0.721, 0.932, 0.994),
  IDis = c(0.898, 0.620, 0.466, 0.621, 0.451, 0.755, 0.318),
  IEve= c(0.46, 0.713, 0.789, 0.644, 0.58, 0.633, 0.543), 
  IUni =c(0.566,0.884, 0.833, 0.599, 0.73, 0.713, 0.758))%>%
  mutate(across(group, factor, levels=c("25m", "370m", "555m", "715m", "1000m", "1335m", "1010m-bottom proximity")))

# plot 
ggradar::ggradar(trophic_indices,
                 values.radar = c("0", "0.5", "1"),
                 group.colours ="#004291", group.line.width = 1,
                 axis.label.size= 4, group.point.size = 3,
                 grid.label.size = 4)+
  facet_wrap(~group)+
  theme_minimal()+
  theme(strip.text.x = element_text(size=12, face = "bold"),
        plot.background = element_rect(colour = "white"))+
  guides(col="none")

Show the code

#ggsave("spider_web.png", path = "figures", dpi = 700)

# spider web for each station 
# for (i in 1:7) {
#  print(ggradar::ggradar(trophic_indices[i,],
#                  values.radar = c("0", "0.5", "1"),
#                  group.colours ="#0c3d97",  group.line.width = 1.3,
#                  axis.label.size= 4, group.point.size = 5,
#                  grid.label.size = 4)+
#   facet_wrap(~group)+
#   theme_light()+
#   theme(strip.text.x = element_text(size=12, face = "bold"))+
#   guides(col="none"))
# }

Show the code

htmltools::tagList(DT::datatable(trophic_indices))

6.additional analyses

Isotopic niches by sations

Show the code

# colors selection 
nichecol <- c("#E4A33A", "#F67451", "#D664BE", "#3DA5D9", "#94b3ae",
              "#18206F", "#FD151B", "#049A8F", "#072AC8", "purple", 
              "#d193f7", "#d8c2ab", "#678FCB", "#A63A49", "#00547A", "#6B54A0")

# plot
ggplot(data = isotope_data_fish, 
                         aes(x = d13c, 
                             y = d15n)) + 
  facet_wrap(~factor(trawling_depth))+
  geom_point(aes(color = species), alpha=0.7, size=1) +
  scale_color_manual(values=nichecol)+  
  scale_fill_manual(values=nichecol)+
  scale_x_continuous(expression({delta}^13*C~'\u2030')) +
  scale_y_continuous(expression({delta}^15*N~'\u2030'))+
  stat_ellipse(aes(group = species, fill = species, color = species), 
               alpha = 0.2, level = p.ell,linewidth = 0.7, type = "norm", geom = "polygon")+
  theme_light()+
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10),
        axis.title = element_text(size=10),
        axis.text = element_text(size=10))+
  labs(shape="Taxon", col= "Species", fill="Species")

Show the code

#ggsave("niches_stations.png", path = "figures", dpi = 700, height = 8, width = 10)

Krill data

• Total \(\delta\)¹⁵N variability = 0.76‰
• Total \(\delta\)¹³C variability = 1.1‰

Show the code

# Load data
isotope_data <-  utils::read.csv(here::here("data", "isotopic_data_2021.csv"), sep = ";", header = T,dec = ",")
isotope_data_krill <- isotope_data %>% filter (species == "Meganyctiphanes_norvegica")

first_plot <- ggplot(data = isotope_data_krill, aes(x = d13c, y = d15n)) + 
  geom_point(aes(color= factor(trawling_depth)), size = 3) +
  ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
  xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
  theme(text = element_text(size=16)) + 
  theme_light()+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")

# Summarise By depth 
data_krill_sum <- isotope_data_krill %>% 
  group_by(trawling_depth) %>% 
  summarise(count = n(),
            mC = mean(d13c), 
            sdC = sd(d13c), 
            mN = mean(d15n), 
            sdN = sd(d15n))

second_plot <- first_plot +
  geom_errorbar(data = data_krill_sum,
                mapping = aes(x = mC, y = mN, ymin = mN - 1.96*sdN, ymax = mN + 1.96*sdN, col = factor(trawling_depth)), 
                width = 0, size=0.8) +
  geom_errorbarh(data = data_krill_sum, 
                 mapping = aes(x = mC, y = mN,xmin = mC - 1.96*sdC, xmax = mC + 1.96*sdC, col = factor(trawling_depth)),
                 height = 0, size=0.8) + 
  geom_point(data = data_krill_sum, 
             aes(x = mC, y = mN, fill = factor(trawling_depth),col = factor(trawling_depth), shape = factor(trawling_depth)),
             size = 5) +
  scale_shape_manual(values = c(21, 22, 23, 24, 25, 8, 7))+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  labs(shape="Trawling depth (m)", col="Trawling depth (m)", fill= "Trawling depth (m)" )

print(second_plot)

Isotope variability

Show the code

boxplot_d13c <- ggpubr::ggboxplot(isotope_data_krill , x = "trawling_depth", y = "d13c", 
          color = "trawling_depth",fill = "trawling_depth", alpha=0.3,
          ylab = "d13c", xlab = "trawling_depth")+
  ylab(expression(paste(delta^{13}, "C (\u2030)"))) +
  xlab("Trawling depth (m)")+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  guides(fill="none", col="none" ,alpha="none")+
  theme_light()

boxplot_d15n <- ggpubr::ggboxplot(isotope_data_krill , x = "trawling_depth", y = "d15n", 
                  color = "trawling_depth", fill = "trawling_depth",
                  alpha=0.3) +
  ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
  xlab("Trawling depth (m)")+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  guides(fill="none", col="none" ,alpha="none")+
  theme_light()

ggpubr::ggarrange(boxplot_d13c, boxplot_d15n)

Tests

Show the code

pairwise.wilcox.test(isotope_data_krill$d13c, isotope_data_krill$trawling_depth,
                 p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum exact test 

data:  isotope_data_krill$d13c and isotope_data_krill$trawling_depth 

     25    370   555   715   1000  1010 
370  0.359 -     -     -     -     -    
555  0.033 0.122 -     -     -     -    
715  0.042 0.125 0.497 -     -     -    
1000 0.033 0.033 0.303 0.179 -     -    
1010 0.033 0.033 0.541 0.125 0.883 -    
1335 0.107 0.433 0.833 1.000 0.454 0.433

P value adjustment method: BH 

Show the code

pairwise.wilcox.test(isotope_data_krill$d15n, isotope_data_krill$trawling_depth,
                 p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum exact test 

data:  isotope_data_krill$d15n and isotope_data_krill$trawling_depth 

     25    370   555   715   1000  1010 
370  0.025 -     -     -     -     -    
555  0.025 0.143 -     -     -     -    
715  0.025 0.243 0.030 -     -     -    
1000 0.025 0.274 0.585 0.025 -     -    
1010 0.025 0.120 0.629 0.025 0.467 -    
1335 0.025 0.306 0.037 0.750 0.025 0.025

P value adjustment method: BH